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Usa  of  the  finite  eluent  aettiod  for  solution  of  soil  consolidation  probleos  is  hcII  established.  Since 
the  first  application  of  this  oethod  [1,21,  considerable  progress  has  been  aade  in  tha  theoretical  foraulation 
as  aall  as  tha  computational  procedures.  Recent  advances  include  variational  foraulations  admitting  liaited 
saoothness  of  finite  eleaant  bases  [3,41,  experimentation  with  several  different  spatial  interpolation  schemes 
C3~91  and  investigation  of  various  temporal  approximation  methods  [6,9-151.  The  finite  element  method  has  been 
applied  to  soils  exhibiting  secondary  compression  [8,11,131,  nonlinear  soil  behavior  [13,16-201,  and  to  finite 
deformation  [16,18,201.  Developments  in  solution  procedures  include  use  of  Laplace  transforms  [10,11,141, 
automatic  selection  of  the  time-step  size  [131  and  the  use  of  single  variable  formulations  [21,221.  Other 
developments  cover  use  of  boundary  element  methods  [23,241,  allouance  for  infinitely  distant  boundaries, 
material  interfaces,  etc. 

One  particular  aspect  of  the  finite  element  procedures  for  soil  consolidation  has  been  the  difficulty  in 
reproducing  accurately  the  pore-nater  pressures  near  loaded  free-draining  surfaces  immediately  after 
application  of  loads  at  these  surfaces.  This  can  be  troublesome  in  problems  involving  inelastic  soil  behavior. 
Taboo's  [231  foraulation,  not  requiring  the  fluid  pressure  to  satisfy  the  prescribed  boundary  condition  at  the 
free-draining  boundary  gave  good  solution  iasediately  after  loading.  Buchaaier’s  [261  studies  shotted  that 
transition  elements,  using  higher  order  interpolation,  near  loaded  boundaries  had  only  liaited  success. 
Vermeer  [271  suggested  that  the  fluid  pressure  boundary  condition  be  enforced  as  a  'ramp'  condition  to  limit 
the  error. 

Herein  ue  report  nuoerical  performance  of  another  approach  to  development  of  finite  element  procedures  to 
satisfactorily  model  the  pore-pressure  distribution  in  the  vicinity  of  free-draining  loaded  surfaces.  It 
consists  of  using  ’singularity’  elements.  In  the  second  section  ue  state  the  tuo-field  foraulation  of  the 
boundary  value  problem  of  consolidation  of  linear  elastic  soils.  The  third  section  contains  a  discussion  of 
the  problem  of  determining  the  pore  pressures  in  the  vicinity  of  a  free-draining  loaded  surface  and  a 
description  of  ’singularity’  elements.  The  numerical  performance  of  several  variants  of  the  scheme  in 
solution  of  one-dimensional  consolidation  is  discussed  in  the  fourth  section.  The  fifth  section  lists  the 
conclusions  arrived  at  as  a  result  of  this  investigation. 


II.  EQUATIONS  60VERNIN6  LINEAR  ELASTIC  SOIL  CONSOLIDATION 


II. A  The  Differential  Equations 

Assuaing  pore  eater  to  be  incoepressible,  the  equations  of  force  equilibrium  of  eleeentary  voluees  and  of 
■ass  continuity  say  be  uritten  in  standard  indicial  notation  as  11,7,281: 

IE  u  1  +7T  +Pf  =0  (1) 

klij  k,l  ,i  ,j  '  j 

IK  OT  ♦  Pi  )1  *u  =0  (2) 

ij  , j  1 2  i  ,i  i,i 
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Cult  )),  Cult  )}  *  victors  of  nodal  point  values  of  the  coaponents  of  displaceeent  at  tiee  t  .  t 
10  10 
respectively; 

dT(t  )),  (HU  )}  s  vKtors  of  nodal  point  values  of  the  pore  water  pressure  at  tiees  t  ,  t 
10  10 
respectively; 

(M  *  the  vector  of  nodal  paint  loads  including  applied  nodal  loads,  boundary  tractions,  body 

forces,  initial  stresses  and  effect  of  displaceeent  constraints; 

(F  )  *  the  vector  of  nodal  point  fluxes  including  applied  nodal  fluxes,  boundary  fluxes,  body 
2 

force  effects  and  effects  of  speciaafied  pore  water  pressures; 

(Cl  *  the  spatial  'stiffness  satrix’  for  the  elastic  soil; 

(Kl  *  the  spatial  'flow  eatrix’  for«<4t  =  1; 

04  *  the  coefficient  characterizing  single  step  teeporal  discretization;  or  ^5  for  stabilit 
(Cl  *  the  coupling  eatrix  representing  the  influence  of  pore  water  pressures  in  the  force 
equilibriui  equation; 

T 

(Cl  *  the  coupling  utrix  representing  the  influence  of  soil  voluee  change  upon  the  nodal  point 
flu. 

The  eatrices  (Kl  and  (Kl  depend  upon  the  interpolation  scheaes  for  displaceaents  and  pore  water  pressures, 
respectively.  The  coupling  eatrix  (Cl  involves  spatial  interpolation  for  both  the  field  variables.  The 
teeporal  discretization  for  the  single  step  scheae  is  reflected  in  the  value  of  the  coefficient 

Equation  (7)  includes  the  'natural’  boundary  conditions  expressed  by  Equations  (4)  and  (6).  Equations  (3) 

^  A 

aed  (3)  are  satisfied  by  uplicitly  requiring  u  *  u  on  S  andir*TTon  S  . 

i  i  li  3 

IL1  fciiiil  MkmUUb 

The  basic  interpolation  scheae  used  in  the  present  study  was  the  8-4  quadrilateral  eleacnt  introduced  by 
Sandhu  [7,301.  In  this  eleaent,  the  displaceaents  have  biquadratic  lagrange  interpolation  and  the  pore  water 
pressure  has  bilinear  lagrange  interpolation  which  is  isoparaaetric  with  the  eleaent  geonetry.  In  this  report, 
this  procedure  is  referred  to  as  the  standard  PS84  procedure.The  singularity  eleaents  used  near  the  loaded 
free-draining  surfaces  will  be  described  in  the  next  section. 
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III.  PORE-PRESSURES  NEAR  LOADED  FREE-DRAININ6  BOUNDARIES 

III. A  Preliainaries 

Iaaediately  after  application  of  a  surface  load  to  a  free-draining  boundary,  the  excess  pore  eater 
pressure  retains  zero  at  the  surface  but  has  a  very  steep  gradient  and  rises,  over  an  extreaely  short  distance 
into  the  soil  aass,  to  a  aagnitude  comparable  with  the  applied  stress.  Finite  eleaent  interpolations  coaaonly 
used  do  not  have  1291  the  capability  to  aodel  this  locally  high  pore  pressure  gradient  near  the  surface. 
Atteapts  to  overcoae  this  difficulty  have  included  use  of  variational  principles  ahere  the  pressure  field  is 
not  required  to  satisfy  the  specified  boundary  condition  1251,  using  fine  aesh  near  the  loaded  surface,  and 
using  higher  order  interpolation  near  the  loaded  surface  1261.  None  of  these  gives  satisfactory  results. 
During  consolidation  of  nonlinear  soils,  the  stress  solution  at  any  tiae  step  defines  the  aechanical  behavior 
for  the  next  step.  Therefoe,  it  is  iaportant  to  get  sufficiently  accurate  eleaent  stresses  at  each  tiae  step. 
In  the  aork  reported  herein,  an  alternative  approach  based  on  use  of  special  'singularity  eleaents’  aas  used. 

III.B  Singularity  Eleaents 

Singularity  eleaents,  using  special  interpolation  scheaes  which  reflect  the  actual  variation  in  the 
variables,  have  been  extensively  used  for  analysis  of  fracture  where  the  stresses  are  unbounded  at  the  crack- 
tip.  Hughes  and  Akin  1311  proposed  special  functions  for  point  as  well  as  line  singularities.  In  the  present 
work  only  line  singularity  was  considered  in  the  context  of  a  two  -diaensional  problea. 

Consider  the  sequence  of  functions 

f  (x)  =  1-  x"  (81 

n 

over  the  doaain  10,11.  Figure  1  shows  the  plots  for  n=0  through  5.  In  the  liait  as  n— f  (x)  =1  for  x  in 

n 

the  interval  10,1)  and  =0  at  x=l.  This  is  the  type  of  discontinuity  encountered  in  consolidation  analysis.  Two 
types  of  singularity  eleaents  were  considered. 

III.Bti  The  two-node  singularity  eleaent 

In  one  diacnsion,  over  range  10,11,  the  singularity  eleaent  would  use  interpolation  functions 

n  n 
1-x  ,  x 

where  n  is  sufficiently  large.  For  n=l  this  reduces  to  linear  interpolation.  Noting  that  the  error  in  pore 
water  pressures  near  loaded  free-draining  boundaries  is  of  relevance  only  iaaediately  after  loading  and  decays 
with  advance  in  the  tiae  doaain,  it  appeared  reasonable  to  use  a  coaposite  eleaent  which  would  approxiaate  the 
singularity  for  seall  values  of  the  'elapsed  tiae’  after  loading  and  reduce  to  linear  interpolation  for  large 
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values  of  the  tiee  variable.  This  led  to  use  of  variants  of  the  type 

n 

f  (x)  *  1  -  ax  -  (1  -  a)x  (9) 

n 

■here  a=0  approximates  the  singularity  eleeent  and  i=i  gives  linear  interpolation.  Thus,  the  coefficient  a  has 
to  be  assigned  a  value  proving  with  tiee  froa  0  to  1.  The  scheee  investigated  ms  based  on 

a  =  1  -  expi-eT)  (10) 

■here  •  is  a  scalar  coefficient  and  T the  non-diaensional  'tiee  factor'.  The  investigation  covered  a  range  of 
values  of  n  and  a. 
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V.  V 


lll.B.ii  The  three-node  sinoularitv  eleeent 
This  eleaent  Mould  directly  involve  three  functions,  viz., 
n 

1,  x,  x 

■here  n  is  sufficiently  large.  This  ms  expected  to  include  linear  interpolation  and  approxiaation  of  the 
singularity  at  the  saae  tiae. 

IIl.B.iii  Interpolation  functions  for  the  quadrilateral  singularity  eleaents 
Special  8-4  and  8-6  eleaents  Mere  used  to  aodel  the  singularity.  Figure  2  shorn  the  arrangeaent  for  the 
tao  eleaents.  The  line  singularity  iaplies  singularity  in  one  variable  only.  In  the  following  m  assuae  this 
to  be  the  variable  t.  For  the  8-4  eleaent  the  interpolating  functions  are: 


{»  *  (1-sMl-t  ) 


n 

s(l-t  ) 


■here  the  range  of  s,  t  is  10,11.  If  the  variant  expressed  by  Equation  (9)  is  eaployed; 
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IV.  NUMERICAL  PERFORMANCE 


IV. A  Introduction 

Several  investigators  133  thru’  39]  have  used  the  finite  eleaent  nethod  to  obtain  approxiaate  solutions 
to  the  problea  of  soil  consolidation.  Using  different  interpolation  scheaes  all  have  generally  reported 
success  uith  ahatever  scheae  they  used.  Comparative  evaluations  of  different  scheaes  are  rare.  Soae 
comparisons  of  nuaerical  pcrforaance  acre  atteapted  by  Sandhu  17,81.  In  evaluating  various  procedures,  Sandhu 
1291  proposed  that  the  folloaing  criteria  be  used. 

i.  The  interpolation  scheae  aust  conform  aith  the  assumptions  regarding  continuity  and  differentiability 
used  in  setting  up  the  governing  variational  foraulation. 

ii.  It  should  be  possible  to  generate  the  'undrained’  solution  i.e.,  the  state  of  fluid  pressures  and 
displacements  at  time  t=0+. 

iii.  For  sufficiently  small  time  steps,  the  scheme  should  be  insensitive  to  the  choice  of  the  time  step 

size. 

The  8-4  element  is  knoan  to  satisfy  all  the  three  conditions  17,8,291  except  that  near  drained  loaded 
surfaces,  the  solution  at  t  3  0*  has  an  error  ahich  decays  aith  advance  in  the  time  domain  171.  If  the 
solution  immediately  after  loading  aas  of  no  interest,  the  8-4  interpolation  scheae  mould  be  adequate. 
Hoaever,  if  accuracy  of  the  solution  at  t*0*is  important  (e.g.  in  nonlinear  or  cyclic  consolidation), 
modifications  to  the  eleaent  interpolation  scheme  are  necessary.  Indeed  this  aas  the  primary  motivation  for 
the  present  investigation. 

IV. 8  Jbf  Example  Problem 

The  procedures  described  mere  applied  to  Terzaghi’s  problea  of  one  dimensional  consolidation.  For  this 
problem,  tha  theoretical  solution  is  knoan  and,  therefore,  precise  comparisons  mere  possible.  The  dimensions 
of  the  consolidating  soil  column  and  the  soil  properties  mere  the  same  as  in  Sandhu’s  example  171.  Figure  3 
shoes  the  soil  column.  Figure  4  shoos  the  finite  eleaent  model.  The  tiae  domain  aas  partitioned  as  tabulated 
boloa: 

1  step  of  A  t  *  0  to  got  the  'undrained’  solution 

1  step  of  At  *  0.00001  over  10.0,0.000011 

1  step  of  A  t  *  0.01  over  10.00001,0.010011 

9  steps  of  At  *0.01  over  10.01001,  0.100011 

10  steps  of  At  *0.1  over  10.10001,  1.100011 

10  steps  of  At  *1.0  over  11.10001,  11.100011 

9  steps  of  At  *10.0  over  111.10001,  101.100011 

10  steps  of  At  *100.0  over  1101.10001,  1101.100011 
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8  steps  of  At  -1000.0  over  [1101.10001,  9101.10001] 


K 


ti  S 


■V 


IV. C  Results  of  the  Analysis 

The  ei aspic  problee  uas  solved  using  the  special  singularity  eleaents  described  earlier.  The  solution  in 
each  case  gave  nodal  point  values  of  the  pore  pressures.  The  error  in  the  solution  ms  non-disensionalized 
through  division  by  the  applied  surface  load. 

Figure  SA  and  Tables  IA,IB  shoe  the  history  of  pore  Mter  pressure  using  the  8-4  eleaent  with  functions 
of  the  type  expressed  by  Equation  (8)  at  depth  .03H  belou  the  top  free-draining  boundary  of  the  soil  coluan  of 
height  H.  Results  stum  that  the  'ioaediate’  pore  Mter  pressure  obtained  by  this  aethod  ms  aore  accurate  than 
the  that  obtained  using  the  PS84.  Houever,  the  pore  pressure  dropped  sharply  after  application  of  the  load  and 
accuracy  in  pore  pressures  ms  uorse  than  that  of  PS84  in  later  tiee  stages.  Apparently,  the  error  in  pore 
pressure  at  later  tiae  stages  uas  due  to  the  absence  of  the  linear  tera  froa  the  interpolation  scheae.  For 
tiae  factor  values  upto  .03  the  error  reduced  uith  n  increasing  upto  3.  For  tiae  factor  values  greater  than  .3 
the  error  increased  uith  increasing  n. 

Figure  SB  and  Table  IC  shoe  the  results  at  a  point  .06H  belou  the  loaded  surface.  Figure  SC  and  Table  IB 
give  the  results  at  a  depth  of  .09H.  As  uould  be  expected,  the  error  at  these  points  is  throughout  SMller 
than  at  .03H.  The  pattern  of  dependence  upon  n  and  the  value  of  the  tine  factor  is  essentially  .pn  23  the 


To  coabine  the  ability  of  the  singularity  interpolation  to  give  accurate  'iaaediate'  pore  pressures  and 
the  accuracy  of  the  PS84  eleaent  for  later  tiae  stages,  it  appeared  reasonable  to  set  up  shape  functions  such 
that  the  eleaent  has  the  character! sties  of  the  singularity  eleaent  at  early  stages  and  these  change  to  those 
of  PS84  as  tiae  increases.  For  this  reason,  shape  functions  of  the  type  expressed  by  Equation  (?)  sere 
developed.  Figures  6A  thru’  6F  and  Tables  2A  thru’  2C  shoo  the  effect  of  variation  in  a  for  values  of  n  equal 
to  3,10,15,30,40  and  48.  Figures  66,  6H  and  Tables  20,  2E  give  additional  inforaation,  for  n  equal  to  48,  for 
the  locations  .06H  and  .09H  belou  the  svface. 

Figure  7  coapares  the  results  for  various  values  of  n  paired  uith  the  optiaaKout  of  the  set  tried)  value 
of  a  in  each  case.  The  least  error  in  initial  porn  pressures  ms  realized  for  n«48,  o*40.  Hoaever,  the  error 
uas  seen  to  grou  uith  tiae.  The  best  overall  accuracy  ms  obtained  for  n«30,  e*32.3.  For  this  case,  the 
aaxiaiua  error  ms  less  than  three  percent.  Figures  8A  thru’  80  shoo  the  spatial  variation  of  error  in  the 
calculated  pore  pressures  for  this  choice  of  n,a. 

Figures  9A  thru’  90  and  Tables  3A  thru’  30  illustrate  the  results  obtained  using  the  8-6  singularity 
eleaent.  Accuracy  of  the  pore  pressure  at  any  location  of  the  entire  soil  coluan  except  .03H  and  .013H  belou 
the  loaded  surface  is  good.  At  later  tiae  stages  also  the  pore  pressures  coapared  sell  uith  the  exact  solution 
and  did  not  differ  each  froa  the  results  of  the  P884  eleaent  regardless  of  the  value  of  n  (so  long  as  it  ms 
greater  than  31  in  the  singularity  eleaent;  n«10  gave  the  best  results. 
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Exact  evaluation  of  the  (learnt  aatricei  requires  the  Gaussian  quadrature  to  use  a  nuaber  of  integration 
points  equal  to  the  index  n  (for  n  an  integer).  This  can  be  quite  expensive.  The  CPU  tiee  using  <8  integration 
points  for  only  one  eleaent  aas  about  50  percent  aore  than  that  for  the  PS84  systea.  To  reduce  coaputational 
costs,  use  of  reduced  integration  order  for  the  8-4  eleaent,  using  interpolants  in  Equation  (8)  aas 
investigated.  The  error  in  pore  pressures,  using  five  point  integration,  increased  with  n  upto  24.  Beyond  this 
value  of  n  the  error  decreased.  The  best  results  Mere  obtained  for  n  equal  to  42.5.  The  results  are  suaaarized 
in  Figures  10A  thru’  10C  and  Tables  40  thru’  40. 

V.  CONCLUSIONS 

Nuaerical  perforaance  of  the  several  scheaes  that  Mere  investigated  indicates  the  folloMing: 

n 

i.  The  8-4  singularity  eleaent  using  shape  functions  based  f(x)  *  l-x  gives  good  values  for  the 

’iaaediate’  pore  pressures  but  at  later  tiae  stages  the  results  are  in  error.  Accuracy  in  pore  uater  pressures 

at  all  tiae  stages  aas  i (proved  using  reduced  order  of  quadrature  for  sufficiently  high  value  of  n. 

n 

ii.  The  8-6  'singularity'  eleaent  using  shape  functions  based  on  f(x)=a*bx+cx  ,  yielded  good  values  of 
the  pare  pressures  at  all  tiae  stages  for  aost  of  the  soil  coluan.  Houever,  at  locations  close  to  the  loaded 
surface,  the  error  aas  quite  large.  This  aakes  the  scheae  unsatisfactory. 

iii.  Accuracy  in  pore  pressures  at  all  tiae  stages  can  be  achieved  using  the  8-4  ‘singularity’  eleaent 

n 

■hose  shape  functions  are  based  an  f (x)  *  1  -  ax  -  (l-a)x  ,  in  ahich  the  coefficient  a  *  1  -  exp (-a  f ) .  This 

ciMbines  tbe  best  characteristics  of  the  singularity  eleaent  Mith  those  of  the  PS84  eleaent. 

The  foregoing  conclusions  are  based  on  study  of  a  one-dieensional  problen.  It  is  inport ant  that  the 

procedures  be  extended  to  probleas  of  tao  and  three  diaensions  (point  singularities  as  aall  as  surface 

singularities)  before  fira  recoaeendations  for  routine  use  of  certain  eleaents  can  be  aade.  (tore  investigation 

is  also  needed  in  tbe  selection  of  indices  n  and  a  used  in  the  foraulation.  Further,  realizing  that  the 

consolidation  of  soils  is  a  decay  process,  use  of  aore  that  one  exponential  teras  in  the  tiae  doaain  could 

possibly  enhance  accuracy.  Use  of  reduced  integration  in  conjunction  aith  interpolation  of  the  type  f(x)  *  1  - 
n 

ax  -  (1  -  als  ,  if  successful,  Mould  reduce  the  coaputational  costs. 
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FIGURE  5A.  DISTRIBUTION  OF  THE  ’RELATIVE'  ERROR  IN 
PORE  PRESSURE  (Special  8-4  Element  , 
Shape  Function  based  on  f(x)*  l-xn  ) 
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FIGURE  7  DISTRIBUTION  OF  THE  ’RELATIVE*  ERROR  IN 
PORE  PRESSURE  (Special  8-4  Element  , 

Shape  Function  based  on  f (x)»l-ax-(l-a)xn  , 


FIGURE  8R  SPRTIRL  DISTRIBUTION  OF  THE  ’RELRTIUE  '  ERROR  IN  PORE  PRESSURE 
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FIGURE  9A  DISTRIBUTION  OF  THE  ’RELATIVE’  ERROR  IN 
PORE  PRESSURE  (Special  8-6  Element  , 
Shape  Punction  based  on  f (x) =a+bx+cxn) 
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FIGURE  10A  DISTRIBUTION  OF  THE  ’RELATIVE’  ERROR  IN 
PORE  PRESSURE  (Special  8-4  Element  , 
Shape  Function  based  on  f(x)»  l-xn, 
5-point  integration  ) 
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FIGURE  10B  DISTRIBUTION  OF  THE  ’RELATIVE *  ERROR  IN 
FORE  PRESSURE  (Special  8-4  Element  , 
Shape  Punction  based  on  f(x)»  l-xn  , 
5-point  integration) 
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